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ABSTRACT 

The onset of runaway stellar collisions in young star clusters is more likely to ini- 
tiate with an encounter between a binary and a third star than between two single 
stars. Using the initial conditions of such three-star encounters from direct A^-body 
simulations, we model the resulting interaction by means of Smoothed Particle Hy- 
drodynamics (SPH). Our code implements new equations of motion that allow for 
efficient use of non-equal mass particles and is capable of evolving contact binaries for 
thousands of orbits. We find that, in the majority of the cases considered, all three 
stars merge together. In addition, we compare our SPH calculations against those 
of the sticky-sphere approximation. If one is not concerned with mass loss, then the 
sticky sphere approach gives the correct qualitative outcome in approximately 75% of 
the cases considered. Among those cases in which the sticky-sphere algorithm identi- 
fies only two particular stars to collide, the hydrodynamic calculations find the same 
qualitative outcome in about half of the instances. If the sticky-sphere approach deter- 
mines that all three stars merge, then the hydrodynamic simulations invariable agree. 
However, in such three star mergers, the hydrodynamic simulations reveal that: (1) 
mass lost as eject a can be a considerable fraction of the total mass in the system (up to 
~ 25%); (2) due to asymmetric mass loss, the collision product can sometimes recieve 
a kick velocity that exceeds 10 km/s, large enough to allow the collision product to 
escape the core of the cluster; and (3) the energy of the ejected matter can be large 
enough (up to - 3 X 10^° erg) to remove or disturb the inter cluster gas appreciably. 



1 INTRODUCTION 



Stars are born in clusters, which upon formation are gener- 
ally dense and massive. In recent years it has become clear 
that clusters remain bound even after losing a considerable 
fraction of their mass due to primordial out-gassing ( Baum- 



gardt Kroupa)2 007 ) . The subsequent dynamical evolution 



collapse evolution, a cluster may enter a phase of gravother- 
mal oscillations ( Cohn et al.|1989 ), allowing periods of high 
interaction rate and providing further opportunity for bina- 
ries and single stars to interact closely. 

Analytic expressions describing encounters between a 
binary and a third star, all treated as point masses, have 



of these clusters leads to a state of core collapse (Portegies 
Zwart et al.|2007 ), almost irrespective of the number of pri- 



been derived for various portions of parameter space ( Heggie 
1975| |Hut|ri983| [Heggie Hut||1993t . In addition, comple- 



mordial binaries (Portegies Zwart et al. 2004); primordial 



binaries do, however, appear to delay the collapse of the 
core ( Fregeau et al.|[2003 Heggie et al.||2006 ). In addition, 
clusters with appropriate initial conditions may form a very 



mentary numerical surveys have been performed in the point 
mass approximation by a number of authors ( [Harrington] 
T970l|Hut k Bahcall|198'3l|Hills|1992] ). During triple encoun- 



massive star by means of runaway stellar collisions (Porte- 
[gies Zwart e t al.|2004 '). Such an object has been hypothised 
to be a progenitor of an intermediate mass black hole (how- 
ever, S( 



Glebbeek et al. (2009)). 



ters, however, individual stars may approach close enough 
to each other that the approximation of point-particle dy- 
namics breaks down: the size and internal structure of the 
stars then play a major role in determining the outcome of 
the encounter. Consequently, some numerical studies have 
augmented the point-mass treatment with simplified mod- 



Even if binaries are not present at the birth of a star 
cluster, they can form via three-body encounters during the 
process of core collapse. Indeed, the expansion of the clus- 



els that incorporate several hydrodynamic effects ( McMillan 



1986; Fregeau etah2004). Large-scale A^-body simulations 
of clusters have demonstrated the ubiquity of resonance in- 



ter core after deep gravothermal collapse ( Sugimoto & Bet- 
twieser|1983 ) is mediated by binaries, regardless of the pres- 
ence or absence of a primordial population. During post-core 



teractions in dynamically unstable triples (Portegies Zwart 



et al. 1999)-the scenarios that ultimately may lead to the 



coalescence of all three stars ( Fregeau et al.|20Q4 ). The accu- 



2 E. Gaburov, J. Lomhardi and S. Portegies Zwart 



rate modelling of the details under which triples merge, and 
whether or not two or all three stars in an encounter partic- 
ipate in the merger, has a profound consequence for the oc- 



currence of collision runaways ( Portegies Zwart McMillan 
2002 Freitag et al.|2006 ) and whether or not such runaways 



can lead to the formation of binaries among intermediate 



mass black holes (Giirkan et al. 2006) 



The first three-dimensional hydrodynamic calculations 
of encounters between a binary and a single star were per- 



formed by Cleary & Monaghan (1990) with the smoothed 



particle hydrodynamics (SPH) method. However, computa- 
tional constraints at that time limited their work to a very 
small number of SPH particles (usually 136 per star) and 
to n = 1.5 poly tropes, appropriate only for white dwarfs 
or extremely low mass main sequence stars. Subsequent hy- 
drodynamic treatments of three-body interactions typically 
confined themselves to scenarios in which at least one of the 
stars was a compact object and therefore could be treated 
as a point mass (e.g.,|Davies et al.|1993l|1994| ). [Davies et al. 



( 1998 ) and Adams et al. (2004) consider three-body encoun- 



ters between a binary and a red giant star as a mechanism 
for destroying red giants near the centres of dense stellar 
systems. Their hydrodynamic simulations follow the fluid of 
the red giant envelope during the encounter, with the red gi- 
ant core and both components of the binary being treated as 
point masses. Because only the red giant envelope is treated 
hydro dynamically, the only mergers that can result are those 
which form a binary of the two point masses surrounded by 
a common envelope donated from the red giant envelope. 
Numerous hydrodynamic simulations of colliding stars 



have studied the structure of the merger product ( 


Benz & 


Hills||1987||Davies et al.||1994| 


Lombardi et al.||1995 


Davies 


et al.||1998| |Lombardi et al.| 


20021 IFreitag & Benz| 2005| 


Gaburov et al.||2008|. In some cases the evolution of these 



collision products is studied further (Suzuki et al. 2007 
Glebbeek et al.||2009 ), especially within the context of the 



formation and evolution of blue stragglers ( Sills et al.|2001 



2005 2008 Glebbeek et al. 2008 Glebbeek k Pols 2008 ) 



Such collision studies, however, have been focused on en- 
counters between two single stars, ignoring for the time be- 
ing that coUisional cross sections and rates can be large for 
systems consisting of three or more stars. 

The scenario of triple-star mergers among low-mass 
main-sequence stars has been previously considered by |Lom-| 
bardi et al. (2003) using SPH. Their calculations indicate 



that the collision product always has a significantly en- 
hanced cross-section and that the distribution of most chem- 
ical elements within the final product is not sensitive to 
many details of the initial conditions. They, however, con- 
centrated solely on low mass stars and treated the triple star 
merger as two separate, consecutive parabolic collisions. 



Recently, Gaburov et al. ( 2008| performed an extensive 
and detailed study to investigate the circumstances under 
which a first collision between stars occurs. Using direct N- 
body integration with sticky spheres of realistic stellar sizes, 
they argued that binaries tend to catalyse collisions. In their 
simulations the binaries that are formed during core collapse 
tend to interact with an incoming star, which subsequently 
merges with one of the binary components. The results of 



Lombardi et al. (2003) suggest that the hydrodynamics of 
such interactions are unlikely to keep the binary itself un- 
damaged. Instead, it is quite likely that the stellar material 



that is expelled during a collision engulfs the system in a 
common envelope, leading to the merger of all three stars. 

In this paper, we introduce a new implementation of 
SPH and apply it to follow accurately the hydrodynamics of 
encounters between hard binaries and intruders. We concen- 
trate on cases involving massive main-sequence stars, such 
as those found in young star clusters, treating all three stars 
simultaneously and with realistic orbital parameters deter- 
mined from a dynamical cluster calculation. In particular, 
the initial conditions are selected from the set of A/'-body 
simulations carried out by Gaburov et al. (2008), but with 



the internal structure of the stars now being determined by 
a stellar evolution code. A comprehensive survey of triple- 
star collisions would need to explore an enormous amount 
of parameter space, but here we focus on a number of rep- 
resentative cases. In total, we selected 40 encounters from 



the simulations of Gaburov et al. (2008). Among these are 
random selections, as well as some that are specifically cho- 
sen because of their relevance for the subsequent A/'-body 
evolution or because of their uncertain outcome given the 
relatively simple treatment of mergers in the A/'-body simu- 
lations. 

This paper is structured as follows. In ^we introduce 
our new formulation of SPH, which allows efficient use of 
non-equal mass particles, as well as our approach for relax- 
ing single stars. In ^we describe how we model close and 
contact binary star systems, and we demonstrate the stabil- 
ity of these systems for at least the time interval of interest. 
The set of initial conditions for the three-body collisions are 
presented in Q Finally, ^presents, while ^discusses, the 
results of our calculations. 



2 METHODS AND CONVENTIONS 
2.1 SPH code 

Smoothed Particle Hydrodynamics is the most widely used 
hydrodynamics scheme in the astrophysics community. It 
is a Lagrangian particle method, meaning that the fluid is 
represented by a finite number of fluid elements or "parti- 
cles." Associated with each particle i are, for example, its 
position r^ , velocity , and mass rrii . Each particle also car- 
ries a purely numerical smoothing length hi that determines 
the local spatial resolution and is used in the calculation of 
fluid properties such as acceleration and density. See [Moii^ 
aghan (1992) and Rasio &: Lombardi| ( |1999| ) for reviews on 
SPH. The code which we used in this work was presented 
in Lombardi et al. (2006). However, we modify the dynami- 
cal equations to allow the efficient use a of range in particle 
masses. This modification is presented in Appendix [A] 



2.2 Choice of units 

Throughout this paper, numerical results are given in units 
where G — Mq — Rq — 1, where G is the Newtonian 
gravitational constant and Mq and Rq are the mass and 
radius of the Sun. The units of time, velocity, and energy 
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Figure 1. Stellar radius versus mass at 2 Myr, as given by the 
TWIN stellar evolution code, for the stars considered in this pa- 
per. 
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2.3 Relaxing a single star 

Before initiating a triple star collision, we must first prepare 
an SPH model for each star in isolation. To compute the 
structure and composition profiles of our parent stars, we use 
the TWIN st ellar evolution code (Eggleton||1971| [Glebbeek 
k Pols|[2008l [Gie bbeek 2 Qd8| from the MUSE software en- 
vironment ( Portegies Zwart et al.|[2009| ^ We evolve main- 
sequence stars with initial helium abundance Y — 0.28 and 
metallicity Z — 0.02 for a time t — 2 Myr, a small enough 
age that even the most massive stars in a star cluster are 
still on the main sequence. The mass-radius relation which 
results from these calculations is shown in Figure [l] 

Initially, we place the SPH particles on a hexagonal 
close packed lattice, with particles extending out to a dis- 
tance only a few smoothing lengths less than the full stel- 
lar radius. After the initial particle parameters have been 
assigned according to the desired profiles from TWIN, we 
allow the SPH fluid to evolve into hydrostatic equilibrium. 
During this calculation, we include the artificial viscosity 
contribution to the SPH acceleration equation so that en- 
ergy is conserved, and we do not include a drag force on the 
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Figure 2. Properties of the SPH model for a 19. IM© star. Pro- 
files are shown as a function of radius, after relaxation for 730 
time units. The frames in the left column show profiles of pres- 
sure P, density p, temperature T (in Kelvin), and mean molecular 
weight ji in units of the proton mass m^, with the dashed curve 
representing results the TWIN evolution code and dots represent- 
ing particle data from our SPH model. The right column provides 
additional SPH particle data: individual SPH particle mass m^, 
smoothing length /i^, number of neighbours A/'jv, and radial com- 
ponent of the hydrodynamic acceleration ahydro (upper data) and 
gravitational acceleration g (lower data). 



particles. For the relaxation calculations of massive stars, 
we do, however, implement a method to keep low mass par- 
ticles from being pushed to large radii: namely, during the 
initial stages of the relaxation, we implement a variation on 
the XSPH method ( Monagha n|1992[ [2002 ) , in which the ve- 
locity used to update positions is the average of the actual 
particle velocity and the desired particle velocity (zero). All 
our relaxed models remain static and stable when left to 
dynamically evolve in isolation. 

This approach allows us to model the desired profiles 
very accurately, and we present an example in Figure [2] 
where we plot desired profiles and SPH particle data for a 
I9.IM0 star. The structure and composition profiles of the 
SPH model closely follow those from TWIN profiles, and the 
model remains stable when evolved dynamically. 



3 PREPARING A BINARY 

In this section, we present our algorithm to model the close 
binary systems that are used in most of the triple star colli- 
sions (see §[5]). The first step in creating a binary is to relax 
each of the two stellar components in isolation, as described 
in the previous section. In the case of detached binaries, we 
place these relaxed stellar models along the x axis with their 
centres of mass separated by the desired separation r. For 
contact binaries, however, we begin with the stars well sep- 
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arated and gradually decrease the semi-major axis until the 
desired separation is achieved, in order to minimise oscilla- 
tions initiated by tidal forces. In all cases, the centre of mass 
of the system remains fixed in space, which we choose to be 
the origin. 

During the binary relaxation process, the positions of 
the particles within each star are adjusted at each timestep 
by simple uniform translations along the binary axis, such 
that the separation between the centres of mass equals the 
desired separation r. Simultaneously, the angular velocity 
Qorb defining the co-rotating frame is continuously updated, 
such that the net centrifugal and gravitational accelerations 
of the two stars cancel exactly: 



t = 8.12 days 



2 V Z^^iTOiXi ^^2"^»^» 



(4) 



where symbolizes a sum over all particles in star j. 

Here, the Cartesian coordinate x is measured along to the 
binary semi-major axis; Vx,i is the acceleration of particle 
i parallel to the axis of the binary in an inert ial frame. A 
centrifugal acceleration is given to all particles such that 
the system approaches a steady state corresponding to a 
synchronised binary. As in the relaxation process of a single 
star, we also include the artificial viscosity contribution to 
the SPH acceleration equation. 

This approach allows us to create close binaries that 
remain in dynamically stable orbits for many hundreds of 
orbits, if not indefinitely. An example is presented in Fig. [S] 
and Fig. |4] In Fig. |3]we plot column densities of a contact 
binary both before and after dynamical evolution through 
over 600 orbits. In Fig. [i] we show time evolution of vari- 
ous energies for the same binary. The epicyclic oscillations, 
with a period of 650 time units, are clearly visible. The fact 
that the epicyclic period is more than an order of magni- 
tude larger than the orbital period of 35 time units under- 
scores how close this binary is to the dynamical stability 
limit. As a binary approaches this limit, the epicyclic period 
would formally approach infinity ( |Rasio k, Shapiro] |1994| ). 
The innermost dynamically stable orbit then marks the tran- 
sition when the squared frequency of the epicyclic oscilla- 
tions passes from a positive to a negative value, so that the 
qualitative behaviour of perturbations changes from oscil- 
latory to exponential. Throughout the calculation, the per- 
turbations remain small and actually damp with time: the 
internal energy U remains constant to within about 0.03%, 
the gravitational energy W to within about 0.008%, and the 
kinetic energy to within about 0.2%. Meanwhile, the total 
energy is conserved within about 0.0004%. 

In another example, we relaxed a contact binary with a 
92.9M0 and a 53.3M0 star with semi-major axis of 43. SR©. 
In Fig [5] we show snapshots for every 50 time units (0.92 
days) of the binary during and after relaxation process. We 
began the relaxation process of a binary with an initial semi- 
major axis of 55.8R0, and we decrease it to 43.8R0in 500 
time units (9.2 days). The top-most left panel shows a bi- 
nary during relaxation at a time of 350 units (6.5 days), and 
the semi-major axis at this time is equal to 48.OR0. It is 
possible to notice commencement of the mass transfer form 
the primary onto the secondary. At the time of 500 units 
(9.22 days), when the semi-major axis becomes 43.8R0, we 
stop the relaxation and dynamically evolve the system in 
the inertial frame. At this time half of the secondary star 
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Figure 3. A contact binary consisting of a 12.2M0 primary and 
a 6.99M0 secondary both at the end of the relaxation (upper 
frame) and after dynamical evolution through more than 600 or- 
bits (lower frame). Colours represent column density, measured 
in g cm~^ on a log scale, along lines of sight perpendicular to the 
orbital plane. 



is already submerged in the fluid of the primary star. By 
the time of 650 units (12 days), the secondary star is com- 
pletely engulfed in the fluid of the primary star. The bottom- 
most right panel shows a binary at the time of 17300 units 
(319 days), and the semi-major axis maintains its value of 
43.8R0. In Fig. |6]we show energy and semi-major axis of 
the binary as a function of time. 
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Figure 5. The relaxation and dynamical evolution of a close binary between 92. QM© and 53. SM© with the semi-major axis equal to 
43.8R0. The orbital period of the binary is 2.78 days (150.6 time units). The calculation switches from a corotating frame to an inertial 
frame at a time of 9.22 days. 



4 INITIAL CONDITIONS 

The parameter space of three-body encounters is immense, 
leaving no hope to be completely covered with SPH simu- 
lations. The approach we take here is to study part of it 
by using the initial conditions obtained from direct A^-body 
simulation. In particular, we take initial condition for three- 



body collisions from Gaburov et al. ( 2008 ) who carried out 



an extensive set of iV-body simulations of young star clus- 
ters. In these simulations the stars were modelled as hard 
spheres with a given mass and corresponding radius. A col- 
lision occurs when two spheres experience physical contact, 
or in other words, when the separation between the cen- 
tre of these spheres is equal to the sum of their radii. This 
treatment of collisions, known as the sticky sphere approxi- 
mation, conserves total mass and momentum. 

In this paper, however, we resolve the stellar structure 
and focus on isolated close three-body interactions. This can 
be justified since usually such interactions last less than a 



year, and therefore local conditions hardly change on such 
a short timescale. All three-body interactions we split in 
two groups: the interaction between a binary and a single 
star, and the interaction between three single stars which 
are in the middle of a resonant interaction. The latter case 
is straightforward to model, as we need to prepare only re- 
laxed single star models, as described in ^ and then assign 
the appropriate initial positions and velocities to each of 
the stars. The actual dynamical interaction process is then 
modelled using the SPH code. 

In the case of an interaction between a binary and a 
single star, we initially relax the binary as described in ^ 
The binary separation is taken from the A'-body simulations. 
Because most of the binaries have separations of a few stel- 
lar radii, tidal circularisation plays an important role, and 
therefore eccentricity of these binaries is nearly equal to zero. 
In some of the cases, the synthetic stellar evolution part of 
A-body calculations predict a binary separation too small 
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Figure 4. Internal energy U, gravitational potential energy W, 
kinetic energy T and total energy E versus time t for the dynam- 
ical evolution of the contact binary shown in Fig. [s] The orbital 
period is 35 time units, while the epicyclic period is 650 time 
units. 




Figure 6. Internal energy gravitational potential energy 
kinetic energy T, total energy E and semi-major axis a versus 
time t for the dynamical evolution of an isolated close binary 
consisting of a 92.9M0 primary and a 53. SM© secondary star. 
All quantities remain within 0.2% of their initial value throughout 
the simulation of more than 100 orbits, highlighting the abilty of 
our code to evolve stably even those binaries in deep contact. The 
small increase in the total energy occurs due to a few low mass 
particles that are escaping to infinity. 




Figure 7. The orientation of a binary and an intruder star 

to be dynamically stable, and in such cases we relax an SPH 
model of the binary near the smallest possible semi-major 
axis such that the binary remains stable or quasi-stable, such 
that the merger time-scale is at least a few thousand time 
units. 

Table [l] lists the initial positions and trajectories in a 
way that is meant to aid the mental visualisation of each 
case: for example, comparing the periastron separation rp,ib 
to the binary semimajor axis ai2 provides an indication 
of where within the binary the intruder strikes. The ra- 
tio ^ib/|^i2| gives a measure of how much energy is being 
brought to the system by the intruder, relative to the binding 
energy of the binary. A negative value of Eih/\Ei2\ implies 
that the intruder star is bound to the binary, otherwise it 
is initially unbound. We note, however, that the magnitude 
of this ratio is much less than unity, which corresponds to 
a nearly parabolic encounter between the intruder and the 
binary star. Indeed, in almost all of these cases the trajec- 
tory of the intruder about the binary is nearly parabolic 
(0.9 < e < 1.1). The last column indicates the angle of ap- 
proach of the intruder toward the binary, with < ^ < 180°. 
More precisely, the angle is the angle between the angular 
momentum vector L12 of the binary and the angular mo- 
mentum vector (L3) of the intruder calculated about the 
center of mass of the binary (Fig. [7|. For example, ^ = 
corresponds to coplanar trajectories with the intruder orbit- 
ing the binary in the same direction (clockwise or counter- 
clockwise) as the binary is orbiting; = 90° corresponds to 
the third star incident on the binary from a direction per- 
pendicular to the orbital plane of the binary; and = 180° 
again corresponds to coplanar trajectories, although now the 
intruder approaches with an angular momentum that is an- 
tiparallel to that of the binary. All of our initial binaries are 
on nearly circular orbits (ei2 < 0.02), with the exception of 
case 249 (ei2 = 0.41). 

We initiate two types of hydrodynamic calculations. 
The first type, which comprises the majority of our calcu- 
lations, consists of a co-rotating binary intruded upon by a 
third star. In these situations, a circular binary is relaxed as 
we described in ^ If it is a contact binary, then the circular 
orbit is maintained, with the orbital plane and phase being 
shifted to match those of the desired initial conditions. If 
the binary is detached, then the velocity of each star is ad- 
justed to give not only the desired orbital orientation and 
phase, but the eccentricity as well. In this way, we account 
for tidal bulging in the binary components. The third star, 
relaxed as described in ^ is initially not rotating and sepa- 
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249 


74.7 


0.11 


0.15 


101 


5.88 


-1.6e-l 


51 


250 


44.0 


31.9 


0.550 


36.4 


3.48 


-2.1e-l 


63 


253 


53.4 


8.55 


0.583 


22.4 


12.5 


-2.5e-2 


35 


256 


33.4 


2.11 


5.84 


18.6 


7.50 


-5.2e-2 


58 


257 


97.3 


24.9 


5.18 


51.5 


2.69 


-2.9e-l 


132 


258 


90.4 


0.55 


0.93 


28.6 


10.2 


-3.5e-3 


86 


259 


55.9 


21.7 


11.4 


26.3 


6.72 


+2.9e-2 


46 


260 


92.9 


53.3 


13.3 


43.9 


54.6 


+5.7e-l 


14 


267 


28.6 


14.2 


19.1 


26.3 


0.28 


-5.3e-l 


126 


298 


56.7 


25.3 


28.1 


26.2 


18.9 


-1.8e-l 


143 


299 


52.3 


16.9 


52.3 


26.2 


0.0 


-0.952e-4 


94 



id 


mi m2 


ms 




[Mo] 





201 


84.1 


0.25 


27.1 


202 


57.9 


0.12 


29.9 


204 


42.4 


11.5 


16.5 


214 


9.49 


16.8 


17.8 


219 


36.6 


9.10 


10.7 


220 


84.3 


68.3 


32.7 


257 


5.18 


24.9 


97.3 


261 


33.9 


13.2 


9.17 


262 


29.3 


31.5 


18.4 



Table 2. The masses, in solar masses, of single stars which par- 
ticipate in the resonance interaction. The first column show the 
case number, while the following columns give the masses of par- 
ticipating stars 



5 RESULTS 

In this section we report on the results of 40 simulations of 
different encounters between three stars. In terms of compu- 
tational time, most of the runs are performed with N ^ 10^ 
and lasted somewhere between one and two weeks on a mod- 



ern P C equipped with an MD-GR APE3 ([Fukushige et al 
T996^ card or an NVIDIA GPU (|Ha mada k Iitaka'" 2007| 
PortegieTZwart et a l.|[2Q07 | IB ellema n et al. 2008; Gab urov| 
et al.|2009| for both self- gravity calculation. Our higher reso- 
lution calculations {N ~ 10^) typically require a few months 
to complete; the total number of integration steps are usu- 
ally between 10^ and 10^. 



Table 1. In the first column, we present the case identification 
number. The second and third columns show the masses of the 
binary components, while the fourth column gives the mass of the 
intruder. The fifth column gives the semimajor axis ai2 of the bi- 
nary. Columns 6 and 7 gives the periastron separation rp^ib and 
eccentricity eib, respectively, of the equivalent two-body Kepler 
orbit between the intruder and the center of mass of the binary. 
Column 8 gives the ratio of the energy Ei\y in this orbit of the in- 
truder and binary to the binding energy \Ei2 \ of the binary itself. 
Column 9 gives the angle ^, in degrees, between the angular mo- 
mentum of the binary and the angular momentum of the intruder 
about the binary (see Fig.|7|. 



rated from the binary by many times the radius of the larger 
star, which allows us to neglect its tidal effects in the initial 
configuration. 

The second type of hydrodynamic calculation involves 
the collision of three individual stars (Table [5]). These also 
represent cases in which a binary is disrupted by an intruder. 
In these cases the three stars are caught in a long-lived reso- 
nant interaction that would be too computational expensive 
to follow entirely with the hydrodynamic code. Each of the 
three stars is first relaxed by the means described in ^ 
Their initial positions and velocities in the collision calcu- 
lation represent a snapshot from the point mass dynamical 
calculation in which the stars were widely separated but 
nearing the end of their resonant interaction. 
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id 


method 


outcome 


speed 


h 


Eej 








[km/sj 




r-1 r\48 1 

[10 "erg] 


201 


pm 


(1,2,3) ^ (1,2),3 


0.8, 352 








ss 


(1,2,3) ^ ({1,3},2) 











14118 


(1,2,3) ^ ({1,3},2) 


< 0.1 


< 0.001 


< 0.1 




28296 


(1,2,3) ^ ({1,3},2) 


< 0.1 


< 0.001 


< 0.1 




113046 


(1,2,3) ^ ({1,3},2) 


< 0.1 


< 0.001 


< 0.1 


202 


pm 


(1,2,3) ^ (1,2),3 


0.7, 529 








ss 


(1,2,3) ^ ({1,3},2) 











6138 


(1,2,3) ^ ({1,3},2) 


< 0.1 


< 0.001 


< 0.1 




11466 


(1,2,3) ^ ({1,3},2) 


< 0.1 


< 0.001 


< 0.1 




22380 


(1,2,3) ^ ({1,3},2) 


< 0.1 


< 0.001 


< 0.1 




91956 


(1,2,3) ^ ({1,3},2) 


< 0.1 


< 0.001 


< 0.1 


203 


pm 


(1,2),3 ^ (1,2,3) ^ (1,2),3 


6.24, 480 








ss 


(1,2),3 ^ (1,2,3) ^ ({1,3},2) 











10398 


(1,2),3 ^ (1,2,3) ^ ({1,3},2) 


< 0.1 


< 0.001 


0.13 


204 


pm 


(1,2,3) ^ (1,3), 2 


40.6, 143 








ss 


(1,2,3) ^ ({1,3},2) ^ {2,{1,3}} 











11946 


(1,2,3) ^ ({1,2},3} ^ {{1,2},3} 


7.2 


0.13 


12. 




60024 


(1,2,3) ^ ({1,2},3} ^ {{1,2},3} 


3.8 


0.082 


12. 


206 


pm 


(1,2),3 ^ (1,2,3) ^ (1,2),3 


60.7, 136 








ss 


(1,2),3 ^ (1,{2,3}) ^ {{2,3},1} 











14475 


(1,2),3 ^ (1,2,3) ^ ({1,2},3) ^ {{1,2},3} 


1.2 


0.048 


13. 


207 


pm 


(1,2),3 


3.31, 308 








ss 


(1,2),3 (1,{2,3}) 











9492 


(1,2),3 


3.7, 345 








208 


pm 


(1,3),2 


1.02, 555 








ss 


(1,3),2 ({1,2},3) 











15018 


(1,3),2 ^ ({1,2},3) 


< 0.1 


< 0.001 


0.22 


211 


pm 


(1,3),2 


91.4, 824 








ss 


(1,3),2 ^ {1,3},2 


53.1, 203 








11028 


(1,3),2 ^ {1,3},2 


39, 146 


0.026 


6.7 


212 


pm 


(1,2),3 ^ (1,2,3) ^ (1,3),2 


34.3, 139 








ss 


(1,2),3 ^ ({1,3},2) ^ {{1,3},2} 











22080 


(1,2),3 ^ (1,2,3) ^ ({1,3},2) ^ {{1,3},2} 


3.6 


0.17 


66. 


213 


pm 


(1,2),3 


1.82, 724 








ss 


(1,2),3 ^ ({1,3},2) 











13314 


(1,2),3 ^ ({1,3},2) 


< 0.1 


< 0.001 


< 0.1 




125130 


(1,2),3 ^ ({1,3},2) 


< 0.1 


< 0.001 


1.5 


214 


pm 


(1,2,3) (1,2),3 


95, 348 








ss 


(1,2,3) ^ ({1,2},3) 











11016 


(1,2,3) ^ ({2,1},3) ^ {{2,1},3} 


1.4 


0.14 


3.9 


217 


pm 


(1,2),3 ^ (1,2,3) ^ (1,2),3 


0.7, 713 








ss 


(1,2),3 ^ ({1,3},2) 











17442 


(1,2),3 ^ ({1,3},2) 


< 0.1 


< 0.001 


< 0.1 




139656 


(1,2),3 ^ ({1,3},2) 


< 0.1 


< 0.001 


< 0.1 


219 


pm 


(1,2,3) ^ (1,2),3 


49, 255 








ss 


(1,2,3) ^ ({1,3},2) ^ {{1,3},2} 











14160 


(1,2,3) ^ ({1,2},3) ^ {{1,2},3} 


8.1 


0.038 


33. 


220 


pm 


(1,2,3) ^ (1,2),3 


166, 778 








ss 


(1,2,3) ^ ({1,2},3)^ {{1,2},3} 











20178 


(1,2,3) ^ (1,{2,3}) ^ {{2,3},1} 


14. 


0.062 


130 




46296 


(1,2,3) -> (1,{2,3}) -> {{2,3},1} 


11. 


0.062 


130 


222 


pm 


(1,2),3 


38.3, 246 








ss 


(1,2),3 ^ ({1,3},2) ^ {{1,3},2} 











17076 


(1,2),3 


48, 310 








223 


pm 


(1,3),2 ^ (1,2,3) ^ (1,2),3 


43, 454 








ss 


(1,3),2 ^ ({1,2},3) ^ {{1,2},3} 











12456 


(1,3),2 ^ (1,2,3) ^ ({2,1},3) ^ {{2,1},3} 


12. 


0.25 


14. 


224 


pm 


(1,2),3 


0.858, 323 
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ss 

14472 


(1,2),3 ^ 
(1,2),3 ^ 


({1,3},2) 

(1,2,3) ^ ({1,3},2) ^ 


{{1,3},2} 




0.94 


0.023 


3.4 


227 


pm 

ss 

10008 


(1,3),2 ^ 
(1,3),2 ^ 
(1,3),2 ^ 


(1,2,3) ^ (1,3),2 
({1,2},3) ^ {{1,2},3} 
(1,2,3) ^ ({1,2},3) ^ 


{{1,2},3} 


19, 55 

3.6 


0.027 


2.3 


231 


pm 

ss 


(2,3), 1 
(2,3), 1 


(1,2,3) ^ (1,2),3 
{1,2},3 




71, 72 
1.3, 163 








13554 


(2,3), 1 ^ 


(1,2,3) ^ ({1,3},2) ^ 


{2,{1,3}} 


0.25 


0.036 


1.9 


232 


pm 

ss 

13110 


(2,3), 1 ^ 
(2,3), 1 ^ 
(2,3), 1 ^ 


(1,2,3) ^ (1,2),3 
(1,{2,3}) ^ {{2,3},1} 
({1,3},2) ^ {{1,3},2} 




106, 107 

4.0 


0.067 


20. 


233 


pm 


(2,3), 1 ^ 


(1,3),2 




18.7, 487 








ss 

13020 


(2,3), 1 ^ 
(2,3), 1 


({1,2},3) ^ {{1,2},3} 
(1,2,3) ^ ({1,2},3) ^ 


{{1,2},3} 



3.7 


0.17 


20 


236 


pm 

ss 

12672 


(1,2),3 
(1,2),3 ^ 
(1,2),3 ^ 


(1,2,3) ^ (2,3),1 
(1,{2,3}) ^ {{2,3},1} 
(1,2,3) ^ ({1,2},3) ^ 


{{1,2},3} 


39, 58 

5.8 


0.14 


25. 


241 


pm 

ss 

19956 


(2,3), 1 ^ 
(2,3), 1 
(2,3), 1 


(1,2,3) ^ (1,3),2 
({1,3},2) ^ {{1,3},2} 
(1,2,3) ^ ({1,3},2) ^ 


{{1,3},2} 


68.1, 128 

8.0 


0.086 


32. 


242 


pm 


(1,2),3 






6.5, 856 








ss 

10224 


(1,2),3 ^ 
(1,2),3 ^ 


({1,3},2) ^ {{1,3},2} 
(1,2,3) ^ ({1,3},2) ^ 


{2,{1,3}} 




0.26 


0.016 


2.0 


245 


pm 

ss 

16884 


(2,3), 1 ^ 
(2,3), 1 ^ 
(2,3), 1 ^ 


(1,2,3) ^ (1,2),3 
({1,3},2) ^ {{1,3},2} 
(1,2,3) ^ ({1,2},3) ^ 


{{1,2},3} 


31.2, 23.9 

5.3 


0.027 


26. 


246 


pm 


(1,2),3 






3, 194 








ss 
5232 
10554 


(1,2),3 
(1,2),3 ^ 
(1,2),3 ^ 


(1,{2,3}) ^ {1,{2,3}} 

(1,2,3) ^ ({1,3},2) ^ {2,{1,3}} 

(1,2,3) ^ (1,2),3 ^ {1,2},3 



2.5 

11.8, 691 


0.017 
0.010 


7.6 
3.7 




21204 
42294 
84642 


(1,2),3 ^ 
(1,2),3 ^ 
(1,2),3 


(1,2,3) ^ (1,{2,3}) ^ 
(1,2,3) ^ (1,{2,3}) ^ 
(1,2,3) ^ ({1,2},3) ^ 


{1,{2,3}} 
{1,{2,3}} 
{1,{2,3}} 


1.1 
0.71 
3.4 


0.023 
0.020 
0.096 


6.3 
5.3 
7.9 


249 


pm 

ss 

10374 


(1,3),2 ^ 
(1,3),2 ^ 
(1,3),2 ^ 


(1,2,3) ^ (1,3),2 

({1,2},3) 

({1,2},3) 




0.06, 30 


< 0.1 


< 0.001 


< 0.1 




82812 


(1,3),2 ^ 


({1,2},3) 




< 0.1 


< 0.001 


< 0.1 


250 


pm 

ss 

10254 


(1,2),3 ^ 
(1,2),3 ^ 
(1,2),3 


(1,2,3) ^ (1,2),3 
(1,2,3) ^ ({1,3},2) 
(1,2,3) ^ ({1,3},2) 




2.3, 324 

0.4 


< 0.001 


0.54 


253 


pm 


(1,2),3 






2.0, 212 








ss 

10374 


(1,2),3 ^ 
(1,2),3 ^ 


({1,3},2) 
{1,2},3 






1.0, 176 


0.032 


5.3 


256 


pm 

ss 

10200 


(1,3),2 ^ 
(1,3),2 ^ 
(1,3),2 


(1,2,3) ^ (1,3),2 
({1,2},3) ^ {{1,2},3} 
(1,2,3) ^ ({1,2},3) ^ 


{{1,2},3} 


39, 234 

7 


0.15 


11. 


257 


pm 

ss 

10236 


(1,2,3) ^ 
(1,2,3) ^ 
(1,2,3) ^ 


(1,2),3 

({1,3},2) - {{1,3},2} 
({1,3},2) ^ {{1,3},2} 




11, 257 

2.9 


0.087 


36. 


258 


pm 

ss 

10104 


(1,3),2 
(1,3),2 ^ 
(1,3),2 ^ 


({1,2},3) ^ {{1,2},3} 
(1,2,3) ^ ({1,2},3) ^ 


{{1,2},3} 


0.8, 79 

0.6 


< 0.001 


2.0 


259 


pm 

ss 

10272 


(1,2),3 
(1,2),3 ^ 
(1,2),3 ^ 


({1,3},2) ^ {{1,3},2} 
(1,2,3) ^ ({1,3},2) ^ 


{{1,3},2} 


52, 353 


5 


0.085 


30. 


260 


pm 

ss 

22518 


(1,2),3 ^ 
(1,2),3 ^ 
(1,2),3 ^ 


(1,2,3) ^ (1,2),3 

(1,{2,3}) 

{1,2},3 




31.0, 340 


44, 456 


0.024 


6.6 


261 


pm 


(1,2,3) ^ 


(1,2),3 




39, 198 








ss 

10092 


(1,2,3) ^ 
(1,2,3) ^ 


({1,3},2) 
({1,3},2) 





7.1 


0.020 


22. 
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262 


pm 
ss 

10314 


(1,2,3) ^ (1,2),3 

(1,2,3) ^ (1,{2,3}) 

(1,2,3) ^ (1,2),3 ^ {2,1},3 


69, 226 


68, 223 


0.011 


2.1 




pm 

ss 

14934 


(l,3j,z (l,^j,3 

(1,3),2 ^ ({1,3},2) ^ {2,{1,3}} 

(1,3),2 ^ ({1,2},3) ^ {{1,2},3} 


403, 905 

9.6 


0.063 


67. 


298 


pm 

ss 

13818 


(1,3),2 ^ (1,2,3) ^ (1,2),3 

(1,3),2 ^ ({1,2},3) ^ {{1,2},3} 

(1,3),2 ^ (1,2,3) ^ ({1,3},2) ^ {{1,3},2} 


60, 200 

7 


0.21 


52. 


299 


pm 

ss 

10194 
102540 


(1,3),2 ^ (1,2,3) ^ (2,3), 1 
(1,3),2 ^ ({1,2},3) ^ {{1,2},3} 
(1,3),2 ^ ({1,2}, 3) ^ {{1,2},3} 
(1,3),2 ^ ({1,2}, 3) ^ {{1,2},3} 


128, 797 

6.9 
2.7 


0.15 
0.16 


310 
330 



Table 3: Summary of the 40 simulations for three-star interactions. The 
first column gives the case identification number. The second column 
either gives the number N of SPH particles used to simulate this case, 
or names the treatment as "pm" (point mass) or "ss" (sticky spheres). 
The third column summarizes the interaction that resulted by listing all 
changes in the state of the system. The fourth column lists the projected 
speed (s) at infinity of the resulting object (s), in units of km s~^. The 
fifth column gives the fractional mass loss 200 time units after the final 
change of state, while the sixth column lists the total energy ejected, in 
units of 10^^ erg, at that same time. 
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5.1 Selected cases 

In Table |3] we summarise the outcomes of all collisions from 
Tables [l] and [2] Binaries are represented by (1,2) while res- 
onances are represented by (1,2,3), with the masses satis- 
fying Ml > M2 > Ms. The merger product between stars 
1 and 2, due either to a binary coalescence or a direct col- 
lision, is represented using braces, {1,2}, where the more 
massive component at the time of the merger is listed first. 
In addition, the notation can be embedded. Consider, for 
example, case 242 with the following interaction sequence: 
(1,2), 3 (1,2,3) ({1,3}, 2) {2, {1,3}}. The initial state 
(1,2), 3 represents a primary 1 and a secondary 2 in a binary 
being intruded upon by the least massive star 3. The (1,2,3) 
indicates that there is neither an immediate retreat of the 
intruder nor an immediate merger, but instead the three 
stars move in a resonance interaction. The state ({1,3}, 2) 
means that the intruder has merged with the primary, leav- 
ing the merger product in a binary with the secondary star. 
Finally, {2, {1,3}} indicates that these two remaining objects 
coalesce. Note that in this final state, the secondary star in- 
dicated to the left of {1,3} within the outer braces, because 
the former was more massive at the time of the merger due 
to dynamical mass transfer during the final stages of binary 
inspiral. 

Here, we present several cases in greater details. First, 
in Figures [S] and [9] we show trajectories and column density 
plots respectively from calculations of case 232, in which a 
I2.2 + 6.99M0 binary collides with a 19. IM© intruder. The 
set up of the initial conditions for these three particular 
stars is described in ^ (Figs. [2] [3] and [4|. Figure [9^ shows 
the three bodies shortly after the start of the calculation. 
Figure[9]3 shows the three bodies just prior to the impact and 
merger of the intruder and the secondary from the binary. 
The first apocentre passage in the resulting binary star is 
shown in Figure [9]:, while Figure |9]i shows the binary in the 
process of merger. In Figure|9^ we show the snapshot shortly 
after the fluid from the three stars has merged into a single 
object, and finally. Figure |9f shows a snapshot from near the 
end of our calculation: the merger product has drifted away 
from the origin due to asymmetric mass loss. In this case, 
the merger product has little angular momentum, and the 
calculated mass loss quickly asymptotes to a constant value 
of approximately 2.6Mq (see Fig. 10). 

As another example, we consider case 260 in which a 
massive binary (92.9 and 53.3M0) is perturbed by a less 
massive intruder (I3.3M0). Here, the intruder is a catalyst 
which triggers binary merger. In Figure [ll] we show time 
evolution of energies (left panel) and global quantities (right 
panel), such as the masses of individual stars and ejected 
fluid. In Figure [12] we present time snapshots for this run. 
Figure shows a snapshot at the beginning of the simu- 
lations, and Figure at the moment of closest approach 
between the intruder star and the binary. The binary merger 
process is shown in Figures |12|i and |12^ . It can be seen 
that fluid is gradually lost from the L2 Lagrangian point. Fi- 
nally, the merged binary is shown in Figure p^. In contrast 
to case 232, binary orbital angular momentum is converted 
into spin of the product, explaining the elongated shape of 
the collisions product in Figure [l2f. One may also notice 
that the collision product is quickly drifting away from the 
centre with a velocity of 14 km/s. Most of this kick velocity 
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Figure 8. Trajectories in the xy (lower right), xz (upper right), 
and yz (upper left) planes for case 232, as given by our hydrody- 
namics calculation. The initial conditions are marked by squares, 
while the final position of an object before merger is marked by 
a 5-point asterisk. 




Figure 10. Masses versus time for case 232. The top frame shows 
the time evolution of the mass of the primary, secondary, and in- 
truder. The bottom frame shows the amount of mass ejected (solid 
curve) and bound in a circumbinary envelope (dotted curve). 
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Figure 9. Column density along lines of sight perpendicular to the xy plane at various times for the same hydrodynamic calculation of 
case 232 presented in Figure [s] 



comes form the escaping intruder star rather than from the 
asymmetric mass ejection. 

The second episode of mass ejection, which occurs after 
the binary merger as can be seen in the right panel of Fig- 
ure [TT] is an artifact of the artificial viscosity used in SPH. 
Initially, the collision product is in the state of both differen- 
tial rotation and hydrostatic equilibrium. Artificial viscosity 
tends to transfer angular momentum from the rapidly ro- 
tating shells to slower ones (Lombardi et al.^l999J , and this 
forces the product to be a solid rotator. Since the inner re- 
gions of the collision product are spinning much faster than 
the outer ones, the angular momentum is transferred out- 
wards. The net effect is that the inner regions of the collision 
product contract, because of loss of the rotational support, 
but the outer regions expand because of the continuously 
increasing supply of angular momentum, and this case can 
be seen in Fig. [13] Eventually, these outer regions become 
unbound and escape, and this results in the second episode 
of the mass loss. The mass loss we report in Table ^ is before 
this second episode but after the first, corresponding to the 
plateau Mejecta ~ 4M0 near t — 7500 in Figure [TT] 

From the data of Table [3] it is evident that when only 
two stars merge the mass loss remains below a few percent, 
and often considerably smaller. It is known that mass loss 



in a parabolic collision between two main-sequence stars is 
small dFreitag Benz|2005||Gaburov et al.|2008t . The mass 
loss percentage is typically larger in cases where all three 
stars ultimately merge, exceeding 10% in the hydrodynamic 
simulations of cases 212, 214, 223, 233, 236, 256, and 298. 
The hydrodynamic evolution in these more extreme cases is 
qualitatively similar: the first merger event is between the 
most massive star and one of the other two, and typically oc- 
curs after a short resonant interaction. The resulting merger 
product is enhanced in size by shock heating and rotation, 
leaving its outermost layers loosely bound. The third star, 
often after flung out to a large distance, can experience sev- 
eral periastron passages through the envelope of the first 
merger product before ultimately donating its fluid to the 
mix. In the process, substantial amounts of gas are ejected 
from the diffuse envelope at every periastron passage. 

An example of this type of interaction is summarised 
in Figures [T4] and [T5] for case 298, which involves a 56. TM© 
+ 25.3M0 binary is intruded upon by a 28.1 M© star. The 
features of these curves can be associated with events during 
the encounter. In this situation the intruding star initiates a 
short lived resonance that ends with the induced merger of 
the binary components near t = 400. As can be seen in the 
middle frame of Figure [15] approximately 2Mq of fluid is 
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Figure 11. In the left panel, we show the time evolution of internal energy U, gravitational potential energy W, kinetic energy T and 
total energy E for case 260, while on the right panel we display the evolution of stellar masses and separations. The top frame shows the 
masses of components 1, 2, and 3, represented by blue, green, and red curves, respectively. The middle frame plots the amount of mass 
ejected (solid curve) and bound in a circumbinary envelope (dotted curve). The bottom curves show the separations between components 
1 and 2 (green), between 1 and 3 (red), and between 2 and 3 (black). 



ejected in the process. The intruder retreats on an eccentric 
orbit, reaching an apastron separation of more than 200 Rq 
and returning for its next pericentre passage shortly before 
t = 1000. As the intruder moves through the outer layers 
of the first merger product, its orbit decays and more mass 
is ejected. By t = 1400, the three-body merger product is 
formed and more than 2OM0 has been ejected in total. 

Another double merger resulting in significant mass loss 
is summarised in Figure |16| which shows the masses and 
separations relevant to the hydrodynamic calculation of case 
256 (a 33.4 + 2.IIM0 binary and a 5.84M0 intruder). Here 
the initial merger occurs between the two most massive stars 
near t — 1800, with about I.SMq of fluid is ejected in the 
process. The third star is left on a highly eccentric orbit, 
reaching an apastron separation of more than QOORq at t = 
4100, and returning for its next pericentre passage at t = 
7200. With each passage through the envelope of the first 
merger product, the orbit of the third star decays and more 
mass is ejected until ultimately, at t = 8000, the three-body 
merger product is formed. 

In cases 250 and 261 the impact of the relatively low- 
mass intruder into the primary causes the outer layers of 
the latter to expand and overflow its Roche lobe, resulting 
ultimately in a stable binary. Figure [T7| shows masses and 
separations of stars for case 261, which begins with the three 
stars in a resonant interaction. At t = 83, the lowest mass 
star is absorbed into the largest star. The collision immedi- 
ately ejects IM0 of material and leaves the two remaining 
stars in an eccentric binary (e ~ 0.4). A fraction of a so- 
lar mass is also placed into a circum-binary envelope: this 
fluid is not gravitationally bound to either star individually 
but rather to the remaining binary as a whole. As the binary 
grinds through the envelope, the orbit gradually circularises. 




Figure 14. Internal energy [/, gravitational energy VK, kinetic 
energy T, and total energy E versus time for case 298. Peaks in 
T and associated dips in W correspond to close passes or mergers 
between the stars. Note that the total energy is conserved to 
about 0.2% over the interval shown. 



as can been seen by examining the separation curve in the 

By t ^ 6 X lO'^, the envelope 



bottom frame of Figure 
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6 X 10^ 

has been effectively removed and the binary has essentially 
reached a steady state with an orbital period of 113 time 
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Figure 12. Column density along lines of sight perpendicular to the xy plane at various times for the same hydrodynamic calculation 
of case 260. 



units (50 hours) and a separation of 26Rq. The calcula- 
tion for this case lasted more than 3.3 x 10^ iterations and 
covered a timespan of over 80000 time units (over 4 years 
simulation time). During this calculation, total energy and 
angular momentum were conserved to better than 0.1%. 



5.2 The effect of numerical resolution 

Because of the longevity of three-body interactions, most of 
our simulations are limited to N ^ (1~2) x 10^ particles. 
Even with this relatively low number of particles, a single 
simulation may take a few weeks to complete, as it typically 
needs to span at least several thousand time units. To test 
whether our results are affected by numerical artifacts, we 
recalculated a few of the simulations in high resolution. In 
most cases, the results are only weakly dependent on the 



resolution. In particular, a case of interest is case 204, which 
begins with three single stars in the middle of the resonance 
interaction. In Figure [18] we present the time evolution of the 
energy for two resolutions. One may see from the kinetic en- 
ergy plot that the first close interaction occurs at t ~ 75. 
The further behaviour of the three stars bear characteristics 
of typical resonant interactions, with kinetic and gravita- 
tional potential energy exhibiting aperiodic oscillations of 
different magnitudes until t ^ 200. At this time two of the 
three stars merge (Table [s]) and binary continues to decay. 
In the high resolution case (right panel in Figure 18), the 



merger occurs somewhat earlier than in the low resolution 
case. Because this kind of interaction is chaotic, it is well 
known that the details at the level of trajectories are reso- 
lution sensitive ( [Davies et al. 111993] [Freitag Benz||2005 ). 
However, the final outcome is consistent between the two 
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Figure 13. Column density along lines of sight perpendicular to the xy plane at various times for the same hydrodynamic calculation 
of case 260. While the inner regions of the product become more compact and spherically symmetric, the outer regions increase in size 
and maintain an elongated shape. 
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Figure 15. Masses and separations versus time for the calcu- 
lation displayed in Fig. ^] case 298. The top frame shows the 
masses of components 1, 2, and 3, represented by blue, green, 
and red curves, respectively. The middle frame plots the amount 
of mass ejected (solid curve) and bound in a circumbinary en- 
velope (dotted curve). The bottom curves show the separations 
between components 1 and 2 (green), between 1 and 3 (red), and 
between 2 and 3 (black). 



resolutions: all three stars eventually merge. Moreover the 
mass and energy of the eject a, as well as the kick velocity 
of the merger product, change by at most a factor of two. 
In Figure [19] we show the time evolution of masses of three 
stars, the mass of ejected fluid and the separation between 
stars. 

Another interesting case is 299, in which a massive bi- 
nary (52.3 +16.9 Mq) is intruded upon by a massive star 
(52.3M0). In Figures [2Q| and 



21 



we show the time evolution 
of energies and global quantities, such as the masses of stars, 
the eject a mass, and the stellar separations. Even though 
there are some differences, the general agreement between 



Figure 16. Like Fig. |15| but for case 256. 



these two simulations is excellent. The merger between two 
of the three stars (the intruder and the primary of the bi- 
nary) occurs at t :^ 110, and further binary decay lasts for 
more than 1200 units. Mass loss and energy of ejected fluid 
are consistent between these two runs of different resolution. 

In Figure |22| we examine the effects of resolution for 
four separate simulations of case 202 with the number of 
particles varying by a factor of 15 from the lowest resolu- 
tion treatment to the highest resolution. The agreement is 
excellent, with even the lowest resolution simulation captur- 
ing all important aspects of the orbital dynamics. The small 
bump in the kinetic energy T shortly after the time t — 100 
corresponds to the absorption of the O.I2OM0 star into the 
57.9M0 star, which excites oscillations in the merger prod- 
uct that are visible in the internal energy U and gravitational 
potential energy W curves. The merger product is left or- 
biting the 29.9M0 star in a stable binary with eccentricity 
e = 0.583 and semimajor axis a = 127Rq: the peaks in T 
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Figure 18. Internal energy U, gravitational potential energy W, kinetic energy T, and total energy E versus time t for two simulations 
of case 204 that differ in resolution: = 11946 (left panel) and 60024 (right panel). 
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Figure 19. Masses and separations versus time for two simulations of case 204: = 11946 (left panel) and 60024 (right panel). The 
top frame shows the masses of components 1, 2, and 3, represented by blue, green, and red curves, respectively. The middle frame plots 
the amount of mass ejected (solid curve) and bound in a circumbinary envelope (dotted curve). The bottom curves show the separations 
between components 1 and 2 (green), between 1 and 3 (red), and between 2 and 3 (black). 



and simultaneous dips in W correspond to the periastron 
passages. 

In Figure [23] we show the projected trajectories of the 
three stars in case 246 of masses 42.2, 38.3 and 1.37M0, 
as calculated with a point-mass integrator (top left frame), 
by using sticky spheres (top right frame) and with the hy- 
drodynamics code (bottom four frames) with different res- 
olution. In all cases, the 1.37M0 intruder approaches the 
circular binary on a hyperbolic trajectory with eccentric- 



ity e = 1.09. In the point mass approximation, the intruder 
reaches a minimum separation of 4.9Oi?0 from the secondary 
and then slingshots back outward on a trajectory with ec- 
centricy e — 1.05. The interaction increases the semimajor 
axis of the binary slightly to 28.5i^0, while also perturbing 
its eccentricity to e = 0.0424. In the sticky sphere approx- 
imation, a merger between the intruder and the secondary 
of the binary occurs during the initial pericenter passage. 
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Figure 20. Internal energy U, gravitational potential energy W, kinetic energy T, and total energy E versus time t for two simulations 
of case 299 that differ in resolution: = 10194 (left panel) and 102540 (right panel). 




Figure 21. Masses and separations versus time for two simulations of case 299: = 10194 (left panel) and 102540 (right panel). Line 
types are as in Fig. |19| 



followed shortly thereafter by a second merger with the pri- 
mary. 

The case plays out qualitatively differently when the 
hydrodynamics is followed. The intruder again passes to a 
minimum separation of about 5Rq from the core of the sec- 
ondary, well within its IIRq stellar radius, and then be- 
gins to retreat. The impact, however, transfers energy into 
oscillations of the secondary and the intruder is not mov- 
ing fast enough to escape further than about 40Rq from 



the secondary. The hydrodynamic calculations indicate that 
the intruder makes a second pericenter passage through the 
secondary, but these calculations deviate depending on the 
resolution: the resulting trajectories do not converge as the 
number of particles is increased up to A/" = 84642 due to the 
chaotic nature of the orbits. 

In the case of our relatively low-resolution N — 10554 
calculation of case 246, the intruder is shot out to a distance 
of over IOORq. Finally, the intruder makes one final pass 
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Figure 17. Like Fig. |15|and Fig. |16| but for case 261, and with 
time plotted on a logarithmic scale so that the long term evolution 
and circularization of the resulting binary can be more easily 
observed. 
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Figure 24. Evolution versus time of, from the top of the figure to 
the bottom, the stellar masses, mass in common envelope (dotted 
curve) as well as in ejecta (solid curve), semimajor axis ai2 of 
the binary, and eccentricity ei2 of the binary (t < 4500) as well 
as eccentricity Ci of the third star as it departs from the merger 
product {t > 4500) for the N = 10554 SPH calculation of case 
246. 



through the secondary, and is ejected out of the system on 
a trajectory with e = 1.3. The removal of orbital energy 
from the binary initiates a mass transfer instability. The 
primary canibalizes the secondary and, as the binary merges, 
O.O6M0 of material is ejected. At this time, the blue and 
green curves in Figure [23] merge into a single blue curve 
(see the lower right hand corner of the middle left frame). 
Masses and orbital parameters for this calculation are shown 
in Figure [24] 

The AT = 21204 and 42294 calculations of case 246 yield 
qualitatively similar results. After the third pericenter pas- 
sage of the intruder through the secondary, the two stars 
merge. The resulting binary, surrounded by an envelope of 
gas removed from the secondary by the impacts, ultimately 
merges. In our highest resolution calculation of this case 
(N = 84642), the intruder does not immediately merge with 
either star in the binary, but rather the three stars move 
around one another in a long-lived resonance interaction be- 
fore ultimately all three stars merge. 



Figure 22. Internal energy [/, gravitational potential energy W 
and kinetic energy T versus time t for four simulations of case 
202 that differ in resolution: N = 6138 (bottom curve), 11466 
(second from bottom), 22380 (third from bottom), 91956 (top). 
The energy scale on the left axis corresponds to the low resolution 
N = 6138 case: the other energy curves have been offset by 10, 
20, and 30 energy units to facilitate the comparison. 



6 DISCUSSION AND CONCLUSIONS 

We present a set of hydrodynamical simulations of 40 close 
encounters between three stars. The initial conditions are 
taken from the high-precision direct A'-body simulations of 
Gaburov et al.| ([2008), who studied the onset of collision 
runaway in young star clusters. Most of the collisions (31) 
involve a massive binary star intruded upon by, generally, a 
lower mass star. The rest of the collisions (9) occur between 
three single stars which are in the middle of the resonant 
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Figure 23. Trajectories projected onto the xy plane for case 246 as calculated in a pure point mass approximation (top left), in a 
sticky sphere approximation (top right), by our hydrodynamics code with = 10554 (middle left), with N = 21204 (middle right), with 
= 42294 (bottom left), and with = 84642 (bottom right). We adopt the convention that the trajectory of the most massive star 
is represented by the blue curve, the intermediate mass star by the green curve, and the lowest mass star by the red curve. The initial 
conditions are marked by squares, while the final position of an object before it merges with another one is marked by a 5-point asterisk. 



interaction. All the simulations were carried out with both 
the SPH method and in the sticky sphere approximation. 

If only initial and final states are of interest, the sticky 
sphere method provides the appropriate outcome of the en- 
counter in about 3 out of every 4 cases. In the cases where 
sticky spheres result in a merger between three stars, our hy- 
drodynamic simulations tend to give a similar result. How- 
ever, if one is interested in mass loss, close inspection reveals 
that in a considerable amount of mass can be ejected in dou- 
ble mergers. In addition, the collision product acquires a kick 
velocity, which is usually a result of the asymmetric mass 
ejection. The kick velocity can be sifficietly high to eject the 
merger product to the cluster halo and even to escape. In 
cases where only two stars merge and the third escapes, the 
kick velocity is large enough that the collision product could 
be ejected out of the star cluster completely. Therefore, it is 
not completely unreasonable to expect collision products to 
be observed in the outer regions of young star cluster, and 
the Pistol star in the Quintuplet cluster (Figer et al. 1998) 
may well be a merger product resulting from an encounter 
between a single and a binary star. 

The sticky sphere approximation, however, fails in sev- 



eral cases. On occasion, this approximation predicts the for- 
mation of a binary with a merger product as one of the 
components (cases 214, 253, 260 and 262), an interesting 
outcome from either an observational or theoretical point of 
view. Detailed hydrodynamic modelling of the same cases, 
however, show that a complete merger is a more likely out- 
come, if the interaction is mild; otherwise, the outcome 
is two unbound stars. In another case, the sticky sphere 
method predicts either one (case 207) or two collisions (case 
222) in a system, but the hydrodynamic simulations predict 
a fly-by. These are the cases for grazing encounters which 
result in the ejection of the intruder star. If the semi-major 
axis of the binary is sufficiently large, binaries tend to avoid 
mergers and become eccentric instead. 

For those situations in which the sticky sphere algo- 
rithm predicts a single merger event, the result is incorrect 
in almost half of the situations. It is important to keep in 
mind that the condition for a merger in the sticky sphere 
approximation is energy independent, and therefore if two 
stars with large enough velocities have a grazing collision, 
this method will incorrectly predict a complete merger. 

Thus in an environment with high velocity dispersion. 
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such as galactic nuclei in which the velocity dispersion is 
typically at least an order of magnitude larger than in the 
cores of young massive star clusters, the sticky sphere ap- 
proximation may fail more often. In such environments, the 
merger cross-section is reduced, as grazing interactions be- 



tween stars may not necessarily lead to mergers (Freitag & 
Benz|2005 ) . While this could be improved by a more sophis- 



ticated effective radius of the merger product (we use simply 
Ri -\- R2), it is unlikely that simple recipes can correctly re- 
produce the richness of the hydrodynamic results, especially 
if one is interested in the close interaction between three or 
more stars. 
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All of our collision products posses some amount of an- 
gular momentum. In some cases, the angular momentum is 
large enough that the shape of the collision product sub- 
stantially deviates from spherical symmetry. Evolving such 
an object is a challenge for stellar evolution codes, given 
that even the evolution of non-rotating massive collisions 
product is a formidable task (Glebbeek 2008). In addition, 
there still exist problems on even hydrodynamical grounds, 
as some of our rotating collision products are gradually los- 
ing mass even at the termination of our hydrodynamic cal- 
culations. The reason for this mass loss is due to spurious 
transport of angular momentum outward caused by artifi- 
cial viscosity ( Lo mbardi et al.|1999[ ) , as described in ^ The 
precise timescale of this effect depends on numerical param- 
eters and the treatment of artificial viscosity. For example, 
in case 220, the progression of the stellar collisions is essen- 
tially the same in the N=20178 and N=46296 calculations. 
In the higher resolution simulation, however, the angular 
momentum transport and resulting mass loss in the final 
collision product progresses more slowly. It is worth noting, 
however, that physical angular momentum transport mech- 
anisms, such as stellar winds and magnetic braking, would 
have a similar qualitative effect but on a longer timescale 
( Sills et al.|2005 ). 



Stellar collisions in a young dense star cluster are ex- 
pected to occur in the first few million years of the cluster 



lifetime ( Portegies Zwart et al.||1999 ). At this age, the star 
cluster may still be embedded in a natal gas fLada & Lada 



2003), and therefore if the eject a is energetic enough the 
state of the gas may be considerably disturbed, and such 
mechanism has recently been proposed within the context 
of globular clusters (iUmbreit et al.||2008|). In the case of 
young star clusters, our results suggest that eject a emanat- 
ing from stellar collisions is energetic enough to significantly 
disturb and even eject the remaining gas. Indeed, a young 
massive star cluster with a star formation efficiency of about 
50% has about 10^^ — 10^° ergs in binding energy of the re- 
maining gas. Our results show that the energy of the ejected 
fluid in stellar collisions exceeds 10^^ ergs, and in two cases 
(cases 220 and 299) even 10^° ergs. Since collisions are ex- 
pected to occur in the core of a star cluster, it would be just 
a matter of a few collisions to significantly perturb or largely 
expel the natal gas from the central region. In the case of a 
runaway merger ( Portegies Zwart et al.|[2004 Freitag et al. 
2006 ) , we therefore expect that the gas will be expelled form 



APPENDIX A: DERIVATION OF SPH 
EQUATIONS OF MOTION 

The use of non-equal mass particles in the simulations al- 
lows us to resolve both the core and the envelope of par- 
ent stars. However, during the merger process, particles of 
significantly different mass from two or more parent stars 
mix, and the standard constraint between density and the 
smoothing length, hi = f{pi,Ci), becomes inappropriate. 
Such a constraint naturally involves a constant with dimen- 
sionality of mass, Ci. This constant is usually determined 
during the set up of the initial conditions and therefore re- 
flects the initial mass resolution of particle i, that is the 
initial total mass of the neighbours of that particle. How- 
ever, as the particle migrates from one region to another, 
the mass resolution of the particle should adapt to its new 
environment. If this does not happen, the particle may have 
too few or too many neighbours, depending on whether it 
migrates into a region with, respectively, an average par- 
ticle mass significantly larger or smaller than in its initial 
environment. To mend this, we present a new approach that 
keeps the number of neighbors roughly constant. Here, we 
can draw an analogy with finite- difference hydrodynamics, 
either on fixed or moving meshes: the number of neighbour- 
ing cells that a given grid cell interacts with is also roughly 
constant (exactly constant on a fixed mesh) and is, to some 
degree, independent of the local fluid conditions. 

We propose a continuous constraint between an esti- 
mate of the number of neighbours and the smoothing length. 
Relaxing the condition that the neighbour number estimate 
be an integer, we weight each neighbour with a function that 
depends on its distance from the particle, G{rij/hi), where 
rij is separation between the particle i and the neighbour 
j. Using such a weight function, we estimate the number of 
neighbours of a given particle i as 

Ni = ^G(|r, -r,|,/i,) = ^G„(/i.). (Al) 

j 3 

We find empirically that the following function provides sat- 
isfactory results: 

G{x, h) = V{4h - 4|x - h\,h), (A2) 

where < x < 2/i, otherwise it is equal to zero, and 



Jo 



V{x,h) = 47T / x^W{x,h)dx. 



(A3) 



the central regions before the end of runaway. 



Here, W{x,h) is an SPH smoothing kernel with a compact 
support of 2/i. We use the kernel of [Monaghan Lattanzio] 
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Figure Al. Plotted here versus normalized separation x/h, the 
weighting function G used to help keep neighbor numbers roughly 
constant even when non-equal mass particles are used. 



(|1985| , for which the weighting function G takes on the form 
shown in Figure |A1| Setting Ni to be constant, equation 
(Al) provides the particle with a necessary constraint be- 



tween hi and its estimated instantaneous number of neigh- 
bours in a continuous way. In the calculations presented in 
this paper, we choose Ni = 22, which typically results in 
about 35 to 40 actual neighbors enclosed by the smoothing 
kernel. (It is not surprising that the actual number of neigh- 
bors is consistently larger than the chosen Ni , as can be seen 



by noting from Fig. Al that G < 1.) 

Because our constraint allows particles smoothing 
lengths to be a function of particle coordinates, the varia- 
tional formulation of SPH can be used to derive equations of 
motion ([Springel Hernquist| | 20021 |Monaghan||2002| |Price| 
|fc M onagha n|2007 ). In particular, we consider the SPH La- 



grangian 



Here, mj is the mass of SPH particle j, vj and uj its ve- 
locity and specific internal energy respectively, and c/yj is its 
gravitational potential, which is defined as 

(t)j = ^^mkg{\rj - rk\,hj) = ^^mkgjk{hj), (A5) 



where g{x,h) is the gravitational potential between two SPH 
particles of unit mass. The Euler- Lagrange equations result- 
ing from this Lagrangian are 



' du \ dpj 



1 dSj 



(A6) 



Here, the first term is the hydrodynamic force, m^ah,*, the 
second term is the gravitational force, m^ag^i, and the par- 
tial derivative, {du/dp)s, is evaluated at constant entropy s. 



Using the SPH definition of density, 

Pj ^'^'rnkWdrj - rk\,hj) = "^nikWjkihj 

k k 

we derive its gradient 



(A7) 



dvi 



^^mkViWik{hi)Sij -\- rriiV iWij (hj 



ruk 



dWjk{hj) dhj 
dvi ' 



dhj 



Differentiating Eq. |A1| with respect to we find 
dhj 



Xj 



dvi 



— —^^^^iGjk{hi)6ij — ViGij{hj 



where 



Xj 



dGjk{hj 
dhj 



(A8) 



(A9) 



(AlO) 



With these equations, it is straightforward to derive accel- 
erations due to pressure 



ah 



-E.^.f [^^W^J{h.) - 3^V.ft,(/lO] (All) 
-E,^. J [^^W^J{hJ) - ^V.G.jihj)] , (A12) 

and due to gravity 

ag,2 - y^^j [^i9ij{hi) + Vigijjhj 

3 



Xirrij 

+ ^y^mj ViGij{hj). 
Here, we define two more quantities: 



rrik- 



dhj 



and 



ruk- 



d gik {hi 
dhi 



(A13) 
(A14) 
(A15) 

(A16) 
(A17) 



Following the approach of "MonaghanI ( 2002 ) (see their 



§2.3), we find the rate of change of the specific internal en- 
ergy to be 



dui Pi v-^ , 



ViW^j[hi 



Xirrij 



- \/iGij {hi 



(A18) 

which guarantees conservation of total energy and entropy 
in the absence of shocks. In order to handle shock waves 
while maintaining energy conservation, we augment these 
equations with artificial viscosity ( Monaghan|1997 ) . For the 
calculations of this paper, we implement a variation on the 
artificial viscosity term proposed by Balsara| ( |1995 ): 



\ ri rj / 

with a — 1 and = 2. In our treatment, 

(vi - v^) • (r^ -r^) fi + fj 



Pij — 



r^ - Yj 



Ci -\- Cj 



(A19) 



(A20) 
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if (vi — Vj) • (r^ — Tj) < 0; otherwise fiij — 0. Here Ci is the 



sound speed at particle i. See Lombardi et al. (2006) for the 



definition of the form factor fi and for additional details on 
how the artificial viscosity is encorporated. 

The evolution equations are integrated using a sym- 
plectic integrator with shared symmetrised timesteps, as in 



Springel] ([2005| . Our shared timestep is determined as 

1 



At Min 



in. [(Atij + At^j) 



where for each SPH particle i, we use 

hi 



Ati,i — Cn,i 



Max [Maxj {i^ij) , Maxj 



with 



and 



At2 



1/2 



\dui/dt\ 



(A21) 



(A22) 



(A23) 



(A24) 



For the simulations presented in this paper, Cn,i — 
0.2 to 0.3 and Cn,2 = 0.05. The Maxj function in equa- 



tion (A22) refers to the maximum of the value of its ex- 
pression for all SPH particles j that are neighbors with i. 



The denominator of equation (A22) is an approximate up- 



per limit to the signal propagation speed near particle i. 
The incorporation of At2 enables us to treat shocks without 
drastically decreasing the timestep during intervals in which 
the flow is subsonic. 



APPENDIX B: INITIAL CONDITIONS 

In Table [B] we summarize the raw initial conditions of our 
calculations in order to facilitate comparisons with any fu- 
ture works. 
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Table Bl. The first column gives the case identification number. The second, third, and fourth columns give the masses Mi, M2, and M3 of 
the colliding stars. Columns 5 through 7 and columns 8 through 10 give the position and velocity, respectively, of star 1 in Cartesian coordinates. 
Likewise, Columns 11 through 13 and columns 14 through 16 give the position and velocity of star 2. The position and velocity of star 3 can be 
determined from the constraints that the center of mass be at the origin and that the net momentum is zero. All quantities are in solar units. 
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